import pandas as pd
import numpy as np
import seaborn as sns
import warnings
import pickle
from sklearn import preprocessing
from sklearn.impute import KNNImputer
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import MinMaxScaler
from sklearn.model_selection import cross_val_score
from sklearn.model_selection import RepeatedStratifiedKFold
from sklearn.ensemble import GradientBoostingClassifier
from sklearn.model_selection import train_test_split
from sklearn.model_selection import GridSearchCV
from sklearn import metrics
from sklearn.metrics import *
le = preprocessing.LabelEncoder()
warnings.filterwarnings("ignore")
pd.set_option('display.max_rows',None)
pd.set_option('display.max_columns',None)
Introduction¶
Purpose and general description¶
This work was created with the goal of trying to find a way to predict adverse weather events. Given my lack of domain knowledge, I was asked to treat the problem as "agnostic" as possible: can "math/statistics" alone find a pattern in the data?
Therefore, it becomes necessary to perform this task by answering a precise question, a clear objective, which I identified as follows: can the collected meteorological data be anticipatory of a possible extreme event? At first glance this work lends itself to an unsupervised analysis: however, not having sufficient knowledge to assess the goodness of this alternative, I preferred to bring the problem in the "supervised" side.
To do this I've created a label in correspondence of the days with a quantity of rain greater than 50mm/day, threshold indicated by the Swiss Meteorological Center as "level of attention" for the data area provided (Domodossola).
Since the data are provided with a delay of one day ("today" I get the data of "yesterday") it is not correct and useful to classify those days according to the occurrence of extreme events. Rather I think it is useful to classify the day before the occurrence of the extreme event, to understand if with the data of the previous day we can "anticipate" the occurrence of a particularly adverse event.
Techincal environment¶
This work has been conducted in a "data science" POC fashion, being focused more on the methodologies aspects than on the technical ones. This aspect leaves unresolved many important issues that will have to be addressed in the eventual production phase of the work.
For the time being, it is worth noting that the libraries used for this POC are the standard Sklearn ones, which could make integration with major cloud facilities (AWS, Azure, GCP) easier. In addition, the training phase was built leveraging the Sklearn Pipeline, in order to get as close as possible to the MLOps paradigms (via mlflow, for example). If necessary, further information on versions and libraries actually used can be provided.
Use case¶
This is the use case thought:
1) Data collected from IoT sensor and drones arrives regarding previous day.
2) The classifier is launched to answer the following question: yesterday was one of those days "anticipating" an extreme event?
3) If yes, weather alert activated!
Next steps¶
Future improvements:
1) Use a smaller time range (hours, minutes, etc...).
2) Better pipeline method for capturing detailed aspects (with customize functions).
2.1) Oversampling data with artificial data generator systems.
2.2) Filling empty columns with EM algorithm.
2.3) Rescaling features passing inverse-beta distribution (more useful when data presents highly skewness).
3) Use a time series modeling that can keep track of any seasonality in the data.
4) Implement adversarial validation for checking distribution shift between train and test set.
5) Make it run real time to make prediction from any interval about next 24 hours.
1.0 Data retrival and ETL¶
df = pd.read_excel('0_dataset/DOMODOSSOLA_giornalieri_1988_2021.xlsx')
df.info()
Set the "Date" column as index and correct the formatting of some fields (from object to float).
df = df.set_index('Data')
for i in df:
if df[i].dtypes == 'object':
if i == 'Settore Prevalente':
pass
else:
df[i] = df[i].str.replace(',', '.').astype(float)
Encoding the only categorical variables into a numerical version, saving original NaN's rows.
df['Settore Prevalente'] = df['Settore Prevalente'].astype(str)
df['Settore Prevalente'] = le.fit_transform(df['Settore Prevalente'])
df['Settore Prevalente'] = df['Settore Prevalente'].mask(df['Settore Prevalente']==10,np.nan)
Let's have a look at the pairwise relationship in the data.
sns.pairplot(df)
2.0 Data labelling¶
As written in the introduction of the paper, two new labels were created in order to answer the question posed at the beginning of the paper: are yesterday's data anticipatory of an extreme climate event?
To do this the following steps are necessary:
1) creation of the label in correspondence of the "extreme" event (mean value of rain preciptations between two hours range);
2) creation of the label corresponding to the day before the occurrence of the extreme event (the target of the classifier).
df['precipitazione_media'] = df[['Precipitazione dalle 9 alle 9 ( mm )','Precipitazione dalle 0 alle 0 ( mm )']].mean(axis=1)
df['evento_estremo'] = 0
df['evento_estremo'] = df['evento_estremo'].mask(df['precipitazione_media']>=50,1)
df['anticipatore'] = 0
df['anticipatore'] = df['anticipatore'].mask(df['evento_estremo']==1,1).shift(-1)
df['anticipatore'].value_counts(normalize=False)
Given the large imbalance between classes present, as shown above, it is necessary to treat the problem second:
- anomaly detection, or
- binary classifier after opportune oversampling or undersampling
The use of anomaly detection, being an alternative close to the unsupervised world, I do not consider it for the moment. It's been started with the binary classifier, with the undersampling of the majority class: in this way I avoid to create bias that could occur in the case of oversampling of the minority class (through the creation of data, perhaps, too much artificial).
df_under = df[df['anticipatore']==0]
df_over = df[df['anticipatore']==1]
df_under = df_under.sample(frac=df_over.shape[0]/df_under.shape[0],random_state=0)
df = df_over.append(df_under).sample(frac=1)
3.2 Edit column format¶
df['Settore Prevalente'] = df['Settore Prevalente'].astype(str)
3.2 Data partitiong¶
print('Dimension of the dataset: ',df.shape)
Since we are dealing with a relatively small dataset, it is not advisable to proceed with the creation of a "classic" train, dev and test set.
For this reason will be implemented only the division into train and test set. Then, in order to manage the tradeoff between bias and overfitting, the train set will undergo a stratified cross-validation procedure together with a search for the best hyperparameters using GridSearch.
In order to avoid the occurrence of data-leakage during the cross-validation it is necessary to implement a pipeline that performs the next steps, in an iterative way, for each fold created by the cross-validation procedure. Here are the steps:
1) Filling column holes with NaN via KNN (holes are filled based on the rows closest to the row under consideration)
2) Re-scaling of the variables
3) Classifier (a classic gradint boosting).
X = df.drop(columns=['anticipatore','evento_estremo']).values
y = df['anticipatore'].values
X_train, X_test, y_train, y_test = train_test_split(X, y , test_size = 0.15 , random_state = 1)
Build the Pipeline
steps = list()
steps.append(('imputation',KNNImputer(n_neighbors=5,weights='distance')))
steps.append(('scaler', MinMaxScaler()))
steps.append(('model', GradientBoostingClassifier(random_state=1)))
pipeline = Pipeline(steps=steps)
Parameters of the GridSearch
parameters = {
'imputation__n_neighbors':[3,5,10],
'model__learning_rate':[0.05,0.1,0.5]
}
Training model with stratified cross-validation and GridSearch. The metric choose is the accuracy, since in this first phase the difference between the "cost" of false positives or false negatives is indifferent. Otherwise it would have been necessary to choose between Precision and Recall (or their harmonic mean, known as F1 metric).
k = RepeatedStratifiedKFold(n_splits=3, n_repeats=3, random_state=1)
model = GridSearchCV(pipeline,param_grid=parameters,cv=k,scoring='accuracy').fit(X_train, y_train)
Ouput of the cross-validation
print("Best Estimator: \n{}\n".format(model.best_estimator_))
print("Best Parameters: \n{}\n".format(model.best_params_))
print("Best Test Score: \n{}\n".format(model.best_score_))
print("Mean Test Scores: \n{}\n".format(model.cv_results_['mean_test_score']))
The results of this first test are promising! Before passing the model to the test set it is necessary to verify that there is not the phenomenon of "distribution shift" that occurs when the data destruction of the train set is very different from that of the dev or test set. To verify this phenomenon three alternatives are possible:
1) classical statistical tests (bayesian t-test for example) 2) adversarial validation 3) histogram comparison
For this POC it has been decided to implement the last version as a sort of baseline. If the average result provided by the next function is greater or arounf than 0.5 we can assume that the difference between the distributions of the train and the test set is marginal.
def histogram_intersection(hist_1, hist_2):
minima = np.minimum(hist_1, hist_2)
intersection = np.true_divide(np.sum(minima), np.sum(hist_2))
return intersection
h1 = pd.DataFrame(X_train)
h2 = pd.DataFrame(X_test)
h1 = h1.replace('nan',np.nan)
h2 = h2.replace('nan',np.nan)
h1 = h1.fillna(0.0)
h2 = h2.fillna(0.0)
inters = []
for i in zip(h1,h2):
h1[i[0]] = h1[i[0]].astype(float)
h2[i[1]] = h2[i[1]].astype(float)
max1 = np.max(h1[i[0]])
max2 = np.max(h2[i[1]])
inter = histogram_intersection(h1[i[0]],h2[i[1]])
inters.append(inter)
print('Mean value of the intersection of all the columns is: ',np.mean(inters))
Given the small dataset we're facing, let's assume train and test do not present "too much" distribution shift.
4.2 Testing¶
y_test_predict_grid = model.predict(X_test)
print("Model Test Recall:", metrics.recall_score(y_test, y_test_predict_grid))
print('--------------------------------------------------')
print("Model Test Precision:", metrics.precision_score(y_test, y_test_predict_grid))
print('--------------------------------------------------')
print("Model Test F1:", metrics.f1_score(y_test, y_test_predict_grid))
print('--------------------------------------------------')
print("Model Test Accuracy:", metrics.accuracy_score(y_test, y_test_predict_grid))
print('--------------------------------------------------')
print('Model Confusion Matrix')
cm = confusion_matrix(y_test,y_test_predict_grid,normalize=None)
cmd = ConfusionMatrixDisplay(cm,display_labels=['No Anticipatore','Sì Anticipatore']).plot()
On the test set the positive result of the training phase is confirmed! It is possible to save the model pipeline in order to reuse it every time new data arrive, always in the format present in this POC.
5.0 Output¶
The output of this POC is a trained classifier that can be reused to make further predictions according to the provided data structure.
with open('meteo_adversity_anticipator.pkl','wb') as file:
pickle.dump(model, file)